{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os, pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=2)\n",
    "from collections import Counter\n",
    "import importlib\n",
    "\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.dpi']= 300\n",
    "mpl.rc(\"savefig\", dpi=300)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "\n",
    "%matplotlib inline\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "scriptpath = '..'\n",
    "sys.path.append(os.path.abspath(scriptpath))\n",
    "\n",
    "from cadrres import pp, model, evaluation, utility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dir = '../result/cv_pred/'\n",
    "n_fold = 1\n",
    "indication = 'HNSC'\n",
    "\n",
    "model_spec_name_list = []\n",
    "for m in ['cadrres', 'cadrres-wo-sample-bias', 'cadrres-wo-sample-bias-weight']:\n",
    "    model_spec_name_list += [\"{}_{}\".format(m, indication)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(81, 27)\n"
     ]
    }
   ],
   "source": [
    "gdsc_drug_df = pd.read_csv('../preprocessed_data/GDSC/hn_drug_stat.csv', index_col=0)\n",
    "gdsc_drug_df.index = gdsc_drug_df.index.astype(str)\n",
    "\n",
    "gdsc_drug_list = gdsc_drug_df.index\n",
    "\n",
    "print (gdsc_drug_df.shape)\n",
    "# gdsc_drug_df.head(1).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "drug_log2_max_conc_dict = dict(zip(gdsc_drug_df.index, gdsc_drug_df['log2_max_conc']))\n",
    "drug_log2_median_ic50_dict = dict(zip(gdsc_drug_df.index, gdsc_drug_df['log2_median_ic50']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "42"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdsc_sample_df = pd.read_csv('../data/GDSC/GDSC_tissue_info.csv', index_col=0)\n",
    "gdsc_sample_df.index = gdsc_sample_df.index.astype(str)\n",
    "\n",
    "gdsc_sample_list = pd.read_csv('../data/GDSC/gdsc_all_abs_ic50_bayesian_sigmoid_only9dosages.csv', index_col=0).index.astype(str)\n",
    "\n",
    "indication_sample_list = [u for u in gdsc_sample_df[gdsc_sample_df['TCGA_CLASS']==indication].index if u in gdsc_sample_list]\n",
    "len(indication_sample_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Read predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dict = {}\n",
    "\n",
    "for model_spec_name in model_spec_name_list:\n",
    "    \n",
    "    output_dict[model_spec_name] = []\n",
    "    \n",
    "    for k in range(1, n_fold+1):\n",
    "        pred_dict = pickle.load(open(output_dir + '{}_5f_{}_output_dict.pickle'.format(model_spec_name, k), 'rb'))\n",
    "        output_dict[model_spec_name].append(pred_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Compare obs and pred for all validation sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import f1_score, precision_score, accuracy_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_df_list = []\n",
    "\n",
    "for model_spec_name in model_spec_name_list:\n",
    "    for k in range(n_fold):\n",
    "        \n",
    "        pred_test_df = output_dict[model_spec_name][k]['pred_test_df']\n",
    "        obs_test_df = output_dict[model_spec_name][k]['obs_test_df']\n",
    "        \n",
    "        sample_list_k = sorted(pred_test_df.index[pred_test_df.index.isin(indication_sample_list)])\n",
    "        pred_test_df = pred_test_df.loc[sample_list_k]\n",
    "        obs_test_df = obs_test_df.loc[sample_list_k]\n",
    "        \n",
    "        drug_list = obs_test_df.columns\n",
    "        \n",
    "        results = []\n",
    "        \n",
    "        for i, d in enumerate(drug_list):\n",
    "            \n",
    "            x = obs_test_df[d].values\n",
    "            y = pred_test_df[d].values\n",
    "            sel = ~np.isnan(x)\n",
    "            \n",
    "            x = x[sel]\n",
    "            y = y[sel]\n",
    "            \n",
    "            # spearman\n",
    "            scor, pval = stats.spearmanr(x, y)\n",
    "            \n",
    "            # F1 weighted\n",
    "            x_bool = (x < drug_log2_max_conc_dict[d]).astype(str)\n",
    "            y_bool = (y < drug_log2_max_conc_dict[d]).astype(str)\n",
    "            f1 = f1_score(x_bool, y_bool, average='weighted')\n",
    "            acc = accuracy_score(x_bool, y_bool)\n",
    "            \n",
    "            precent_sensitive = (np.sum(x < drug_log2_max_conc_dict[d]) / len(x)) * 100\n",
    "            \n",
    "            if (precent_sensitive > 0) & (precent_sensitive < 100):\n",
    "                (f1_resistant, f1_sensitive) = f1_score(x_bool, y_bool, average=None)\n",
    "                (precision1_resistant, precision_sensitive) = precision_score(x_bool, y_bool, average=None)\n",
    "            else:\n",
    "                f1_resistant = np.nan\n",
    "                f1_sensitive = np.nan\n",
    "                precision1_resistant = np.nan\n",
    "                precision_sensitive = np.nan\n",
    "            \n",
    "            # MAE (sensitive)\n",
    "            sensitive_sel = x < drug_log2_max_conc_dict[d]\n",
    "            x_sensitive = x[sensitive_sel]\n",
    "            y_sensitive = y[sensitive_sel]\n",
    "            mae = np.mean(np.abs(x_sensitive - y_sensitive))\n",
    "            \n",
    "            results += [[d, precent_sensitive, drug_log2_max_conc_dict[d], drug_log2_median_ic50_dict[d],\n",
    "                         scor, pval, acc, f1, f1_resistant, f1_sensitive, precision1_resistant, precision_sensitive, mae]]\n",
    "            \n",
    "        result_df = pd.DataFrame(results, columns=['drug_id', 'precent_sensitive', 'log2_max_conc', 'log2_median_ic50', 'spearman', 'pval', 'accurary', 'f1', 'f1_resistant', 'f1_sensitive', 'precision1_resistant', 'precision_sensitive', 'MAE'])\n",
    "        result_df.loc[:, 'model'] = model_spec_name\n",
    "        result_df.loc[:, 'k'] = k+1\n",
    "        \n",
    "        result_df_list.append(result_df)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(243, 15)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_result_df = pd.concat(result_df_list, axis=0)\n",
    "all_result_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_result_df.to_excel('../result/cv_pred/cv_score_{}.xlsx'.format(indication), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Summarize 5-fold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>drug_id</th>\n",
       "      <th>precent_sensitive</th>\n",
       "      <th>log2_max_conc</th>\n",
       "      <th>log2_median_ic50</th>\n",
       "      <th>spearman</th>\n",
       "      <th>pval</th>\n",
       "      <th>accurary</th>\n",
       "      <th>f1</th>\n",
       "      <th>f1_resistant</th>\n",
       "      <th>f1_sensitive</th>\n",
       "      <th>precision1_resistant</th>\n",
       "      <th>precision_sensitive</th>\n",
       "      <th>MAE</th>\n",
       "      <th>model</th>\n",
       "      <th>k</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1001</td>\n",
       "      <td>71.428571</td>\n",
       "      <td>10.965784</td>\n",
       "      <td>10.861657</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>0.588724</td>\n",
       "      <td>0.714286</td>\n",
       "      <td>0.714286</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.8</td>\n",
       "      <td>1.095393</td>\n",
       "      <td>cadrres_HNSC</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1003</td>\n",
       "      <td>100.000000</td>\n",
       "      <td>-3.321928</td>\n",
       "      <td>-5.793488</td>\n",
       "      <td>0.600000</td>\n",
       "      <td>0.208000</td>\n",
       "      <td>0.833333</td>\n",
       "      <td>0.909091</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.781029</td>\n",
       "      <td>cadrres_HNSC</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1004</td>\n",
       "      <td>100.000000</td>\n",
       "      <td>-3.321928</td>\n",
       "      <td>-6.119126</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.996345</td>\n",
       "      <td>cadrres_HNSC</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1006</td>\n",
       "      <td>57.142857</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.553154</td>\n",
       "      <td>0.392857</td>\n",
       "      <td>0.383317</td>\n",
       "      <td>0.571429</td>\n",
       "      <td>0.552381</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.6</td>\n",
       "      <td>2.087725</td>\n",
       "      <td>cadrres_HNSC</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1007</td>\n",
       "      <td>85.714286</td>\n",
       "      <td>-6.321928</td>\n",
       "      <td>-7.188990</td>\n",
       "      <td>-0.071429</td>\n",
       "      <td>0.879048</td>\n",
       "      <td>0.571429</td>\n",
       "      <td>0.623377</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.727273</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.8</td>\n",
       "      <td>1.749042</td>\n",
       "      <td>cadrres_HNSC</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  drug_id  precent_sensitive  log2_max_conc  log2_median_ic50  spearman  \\\n",
       "0    1001          71.428571      10.965784         10.861657  0.250000   \n",
       "1    1003         100.000000      -3.321928         -5.793488  0.600000   \n",
       "2    1004         100.000000      -3.321928         -6.119126  0.000000   \n",
       "3    1006          57.142857       1.000000          0.553154  0.392857   \n",
       "4    1007          85.714286      -6.321928         -7.188990 -0.071429   \n",
       "\n",
       "       pval  accurary        f1  f1_resistant  f1_sensitive  \\\n",
       "0  0.588724  0.714286  0.714286           0.5      0.800000   \n",
       "1  0.208000  0.833333  0.909091           NaN           NaN   \n",
       "2  1.000000  1.000000  1.000000           NaN           NaN   \n",
       "3  0.383317  0.571429  0.552381           0.4      0.666667   \n",
       "4  0.879048  0.571429  0.623377           0.0      0.727273   \n",
       "\n",
       "   precision1_resistant  precision_sensitive       MAE         model  k  \n",
       "0                   0.5                  0.8  1.095393  cadrres_HNSC  1  \n",
       "1                   NaN                  NaN  1.781029  cadrres_HNSC  1  \n",
       "2                   NaN                  NaN  0.996345  cadrres_HNSC  1  \n",
       "3                   0.5                  0.6  2.087725  cadrres_HNSC  1  \n",
       "4                   0.0                  0.8  1.749042  cadrres_HNSC  1  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_result_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>drug_id</th>\n",
       "      <th>model</th>\n",
       "      <th>precent_sensitive</th>\n",
       "      <th>log2_max_conc</th>\n",
       "      <th>log2_median_ic50</th>\n",
       "      <th>spearman</th>\n",
       "      <th>accurary</th>\n",
       "      <th>f1</th>\n",
       "      <th>f1_resistant</th>\n",
       "      <th>f1_sensitive</th>\n",
       "      <th>precision1_resistant</th>\n",
       "      <th>precision_sensitive</th>\n",
       "      <th>MAE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1001</td>\n",
       "      <td>cadrres-wo-sample-bias-weight_HNSC</td>\n",
       "      <td>71.428571</td>\n",
       "      <td>10.965784</td>\n",
       "      <td>10.861657</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.857143</td>\n",
       "      <td>0.863492</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.888889</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.026040</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1001</td>\n",
       "      <td>cadrres-wo-sample-bias_HNSC</td>\n",
       "      <td>71.428571</td>\n",
       "      <td>10.965784</td>\n",
       "      <td>10.861657</td>\n",
       "      <td>0.428571</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.002920</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1001</td>\n",
       "      <td>cadrres_HNSC</td>\n",
       "      <td>71.428571</td>\n",
       "      <td>10.965784</td>\n",
       "      <td>10.861657</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>0.714286</td>\n",
       "      <td>0.714286</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.8</td>\n",
       "      <td>1.095393</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1003</td>\n",
       "      <td>cadrres-wo-sample-bias-weight_HNSC</td>\n",
       "      <td>100.000000</td>\n",
       "      <td>-3.321928</td>\n",
       "      <td>-5.793488</td>\n",
       "      <td>-0.085714</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.945527</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1003</td>\n",
       "      <td>cadrres-wo-sample-bias_HNSC</td>\n",
       "      <td>100.000000</td>\n",
       "      <td>-3.321928</td>\n",
       "      <td>-5.793488</td>\n",
       "      <td>0.314286</td>\n",
       "      <td>0.833333</td>\n",
       "      <td>0.909091</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.950149</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  drug_id                               model  precent_sensitive  \\\n",
       "0    1001  cadrres-wo-sample-bias-weight_HNSC          71.428571   \n",
       "1    1001         cadrres-wo-sample-bias_HNSC          71.428571   \n",
       "2    1001                        cadrres_HNSC          71.428571   \n",
       "3    1003  cadrres-wo-sample-bias-weight_HNSC         100.000000   \n",
       "4    1003         cadrres-wo-sample-bias_HNSC         100.000000   \n",
       "\n",
       "   log2_max_conc  log2_median_ic50  spearman  accurary        f1  \\\n",
       "0      10.965784         10.861657  0.500000  0.857143  0.863492   \n",
       "1      10.965784         10.861657  0.428571  1.000000  1.000000   \n",
       "2      10.965784         10.861657  0.250000  0.714286  0.714286   \n",
       "3      -3.321928         -5.793488 -0.085714  1.000000  1.000000   \n",
       "4      -3.321928         -5.793488  0.314286  0.833333  0.909091   \n",
       "\n",
       "   f1_resistant  f1_sensitive  precision1_resistant  precision_sensitive  \\\n",
       "0           0.8      0.888889              0.666667                  1.0   \n",
       "1           1.0      1.000000              1.000000                  1.0   \n",
       "2           0.5      0.800000              0.500000                  0.8   \n",
       "3           NaN           NaN                   NaN                  NaN   \n",
       "4           NaN           NaN                   NaN                  NaN   \n",
       "\n",
       "        MAE  \n",
       "0  1.026040  \n",
       "1  1.002920  \n",
       "2  1.095393  \n",
       "3  1.945527  \n",
       "4  1.950149  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "score_df = all_result_df.groupby(by=['drug_id', 'model']).mean().reset_index().drop(['pval', 'k'], axis=1)\n",
    "score_df.to_excel('../result/cv_pred/cv_score_summary_{}.xlsx'.format(indication), index=False)\n",
    "score_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
